publication: [IL-1β+ macrophages fuel pathogenic inflammation in pancreatic cancer] (https://www-nature-com.proxy.kib.ki.se/articles/s41586-023-06685-2#Sec1)
Published code on [GitHub:] (https://github.com/ostunilab/PDAC_Nature_2023/tree/main/scRNAseq/Mouse/Timecourse_KPC)
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 'SeuratObject' was built under R 4.4.0 but the current version is
## 4.4.1; it is recomended that you reinstall 'SeuratObject' as the ABI
## for R may have changed
## 'SeuratObject' was built with package 'Matrix' 1.7.0 but the current
## version is 1.7.1; it is recomended that you reinstall 'SeuratObject' as
## the ABI for 'Matrix' may have changed
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.4 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(scDblFinder)
## Loading required package: SingleCellExperiment
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
##
## The following object is masked from 'package:dplyr':
##
## count
##
##
## Attaching package: 'MatrixGenerics'
##
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
##
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
##
## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union
##
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
##
## The following object is masked from 'package:SeuratObject':
##
## intersect
##
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
##
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
## tapply, union, unique, unsplit, which.max, which.min
##
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
##
## The following objects are masked from 'package:lubridate':
##
## second, second<-
##
## The following object is masked from 'package:tidyr':
##
## expand
##
## The following objects are masked from 'package:dplyr':
##
## first, rename
##
## The following object is masked from 'package:utils':
##
## findMatches
##
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
##
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
##
## The following object is masked from 'package:lubridate':
##
## %within%
##
## The following object is masked from 'package:purrr':
##
## reduce
##
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
##
## The following object is masked from 'package:sp':
##
## %over%
##
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
##
## Attaching package: 'Biobase'
##
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
##
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
##
##
## Attaching package: 'SummarizedExperiment'
##
## The following object is masked from 'package:Seurat':
##
## Assays
##
## The following object is masked from 'package:SeuratObject':
##
## Assays
set.seed(123)
dirs <- list.dirs(path = './GSE217846_RAW/data/samples/', recursive = F, full.names = F)
for (x in dirs){
cts <- ReadMtx(mtx= paste0('./GSE217846_RAW/data/samples/', x, '/matrix.mtx.gz'),
features = paste0('./GSE217846_RAW/data/samples/', x, '/features.tsv.gz'),
cells=paste0('./GSE217846_RAW/data/samples/', x, '/barcodes.tsv.gz'))
assign(x, CreateSeuratObject(counts = cts, min.cells = 3, project = x))
}
# merge object
seurat_object <- merge(Healthy, y=c(Tumor_d10,Tumor_d20, Tumor_d30),
add.cell.ids=c("Healthy","Tumor_d10", "Tumor_d20", "Tumor_d30"))
# create columns for metadata
# identity
unique(seurat_object@meta.data$orig.ident) # control identity
## [1] "Healthy" "Tumor_d10" "Tumor_d20" "Tumor_d30"
# cell id
seurat_object@meta.data$cell_id <- rownames(seurat_object@meta.data)
# percent mt
seurat_object[["percent.mt"]]<-PercentageFeatureSet(seurat_object, pattern="^mt-")
plot1 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
seurat_object <- subset(seurat_object, subset = percent.mt < 25 & nCount_RNA > 1000 & nFeature_RNA > 200)
head(seurat_object)
## orig.ident nCount_RNA nFeature_RNA
## Healthy_AAACCCAAGTAAACAC-1 Healthy 6472 1961
## Healthy_AAACCCACACTGCTTC-1 Healthy 6203 1950
## Healthy_AAACCCACAGCGACAA-1 Healthy 5495 1804
## Healthy_AAACCCATCGAAACAA-1 Healthy 3126 1095
## Healthy_AAACGAAAGACCTTTG-1 Healthy 1390 234
## Healthy_AAACGAAAGGAGTACC-1 Healthy 4835 1786
## Healthy_AAACGAAAGGAGTATT-1 Healthy 6673 2216
## Healthy_AAACGAAAGGCTAACG-1 Healthy 8308 2464
## Healthy_AAACGAACACTTGAAC-1 Healthy 10632 3141
## Healthy_AAACGAACAGGACAGT-1 Healthy 8371 2322
## cell_id percent.mt
## Healthy_AAACCCAAGTAAACAC-1 Healthy_AAACCCAAGTAAACAC-1 14.292336
## Healthy_AAACCCACACTGCTTC-1 Healthy_AAACCCACACTGCTTC-1 13.557956
## Healthy_AAACCCACAGCGACAA-1 Healthy_AAACCCACAGCGACAA-1 7.770701
## Healthy_AAACCCATCGAAACAA-1 Healthy_AAACCCATCGAAACAA-1 8.925144
## Healthy_AAACGAAAGACCTTTG-1 Healthy_AAACGAAAGACCTTTG-1 10.503597
## Healthy_AAACGAAAGGAGTACC-1 Healthy_AAACGAAAGGAGTACC-1 7.321613
## Healthy_AAACGAAAGGAGTATT-1 Healthy_AAACGAAAGGAGTATT-1 6.069234
## Healthy_AAACGAAAGGCTAACG-1 Healthy_AAACGAAAGGCTAACG-1 7.029369
## Healthy_AAACGAACACTTGAAC-1 Healthy_AAACGAACACTTGAAC-1 6.593303
## Healthy_AAACGAACAGGACAGT-1 Healthy_AAACGAACAGGACAGT-1 9.090909
seurat_object
## An object of class Seurat
## 22211 features across 54311 samples within 1 assay
## Active assay: RNA (22211 features, 0 variable features)
## 4 layers present: counts.Healthy, counts.Tumor_d10, counts.Tumor_d20, counts.Tumor_d30
plot1 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(seurat_object, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
# double check filtering
min(seurat_object@meta.data$nFeature_RNA)
## [1] 201
min(seurat_object@meta.data$nCount_RNA)
## [1] 1001
max(seurat_object@meta.data$percent.mt)
## [1] 24.99137
for (i in c('Healthy','Tumor_d10','Tumor_d20','Tumor_d30')){
sub <- subset(seurat_object, subset = orig.ident == i)
eval(parse(text=paste("sceDblF_",i," <- scDblFinder(GetAssayData(sub, layer = 'counts'), dbr = 0.07)",sep="")))
eval(parse(text=paste("score.",i," <- sceDblF_",i,"@colData@listData[['scDblFinder.score']]",sep="")))
eval(parse(text=paste("names(score.",i,") <- rownames(sceDblF_",i,"@colData)",sep="")))
}
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~6189 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 1267 cells excluded from training.
## iter=1, 1263 cells excluded from training.
## iter=2, 1252 cells excluded from training.
## Threshold found:0.421
## 613 (7.9%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~11676 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2174 cells excluded from training.
## iter=1, 2442 cells excluded from training.
## iter=2, 2470 cells excluded from training.
## Threshold found:0.615
## 1779 (12.2%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~12244 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2233 cells excluded from training.
## iter=1, 2618 cells excluded from training.
## iter=2, 2544 cells excluded from training.
## Threshold found:0.532
## 1770 (11.6%) doublets called
## Assuming the input to be a matrix of counts or expected counts.
## Creating ~13341 artificial doublets...
## Dimensional reduction
## Evaluating kNN...
## Training model...
## iter=0, 2308 cells excluded from training.
## iter=1, 2665 cells excluded from training.
## iter=2, 2745 cells excluded from training.
## Threshold found:0.615
## 1995 (12%) doublets called
doublets.info <- rbind(sceDblF_Tumor_d10@colData,sceDblF_Tumor_d20@colData,sceDblF_Healthy@colData,sceDblF_Tumor_d30@colData)
seurat_object$is.doublet <- doublets.info$scDblFinder.class
seurat_object <- subset(seurat_object, subset = is.doublet == 'singlet')
# default parameters
seurat_object<-NormalizeData(seurat_object)
## Normalizing layer: counts.Healthy
## Normalizing layer: counts.Tumor_d10
## Normalizing layer: counts.Tumor_d20
## Normalizing layer: counts.Tumor_d30
seurat_object<-FindVariableFeatures(seurat_object, selection.method = "vst", nfeatures=2000)
## Finding variable features for layer counts.Healthy
## Finding variable features for layer counts.Tumor_d10
## Finding variable features for layer counts.Tumor_d20
## Finding variable features for layer counts.Tumor_d30
all.genes<-rownames(seurat_object)
seurat_object<-ScaleData(seurat_object, features=all.genes)
## Centering and scaling data matrix
seurat_object <- RunPCA(seurat_object, features = VariableFeatures(object = seurat_object))
## PC_ 1
## Positive: Krt8, Krt18, S100a6, Clu, Lmna, Fbln2, Tm4sf1, Tm4sf4, Krt7, Krt19
## Sox4, Tspan8, Cavin1, Emp1, Cystm1, Onecut2, Nedd4, S100a16, 2200002D01Rik, Cdkn2a
## Cldn18, Anxa2, Scd2, Prkg2, Serpinb6a, Ccnd1, Tpm1, Pmepa1, Msln, Npnt
## Negative: Ctrb1, Try4, Prss2, Cela1, 2210010C04Rik, Cd74, Cela2a, Pnlip, Clps, Reg1
## Tyrobp, H2-Aa, Cpb1, Cela3b, Fcer1g, H2-Ab1, H2-Eb1, Try5, Ctrl, Cpa1
## Zg16, Rnase1, Pnliprp1, Sycn, Bcl2a1b, Il1b, Mpeg1, Ifi27l2a, Irf8, Gm2a
## PC_ 2
## Positive: Tspan8, Krt18, Krt8, Cystm1, Krt7, Onecut2, Tm4sf4, Cldn18, Krt19, Lgals4
## Muc1, Cdh17, Lgals2, Cdkn2a, Msln, Ctse, Cdh1, Prxl2a, Prkg2, Apob
## Npnt, Anxa13, 2200002D01Rik, Cldn6, Krt20, Anxa10, Ccnd1, Dsp, Cldn2, Cpe
## Negative: Aebp1, Fstl1, Serping1, Bgn, Col3a1, Sparc, Col5a2, Loxl1, Dcn, Fbn1
## Col1a2, Meg3, Ccdc80, Col6a2, Serpinf1, Lhfp, Plod2, Adamts2, Slit3, C1s1
## Col1a1, Rcn3, Col6a3, Lox, Cpxm1, Gas1, Serpinh1, Prrx1, Cdh11, Pdgfra
## PC_ 3
## Positive: Il1b, Lyz2, Cxcl2, Fcer1g, Ccl6, C1qa, Ptgs2, Mafb, Tgfbi, C1qc
## C1qb, Fcgr3, Csf1r, Tyrobp, Mpeg1, Thbs1, Apoe, Arg1, Dab2, Msr1
## Ccl9, Clec4n, Ctsd, Pla2g7, Cfp, Fn1, Ms4a6c, C3ar1, Trf, Csf2rb
## Negative: Cela2a, Cpa1, Cela3b, Pnlip, Ctrb1, Try4, Ctrl, Clps, Prss2, 2210010C04Rik
## Rnase1, Sycn, Cela1, Cpb1, Zg16, Reg1, Cel, Try5, Pnliprp1, Cpa2
## Pla2g1b, Uba52, Spink1, Tff2, Gimap4, Rps18, Klk1, Rpl21, Rps15a, Pnliprp2
## PC_ 4
## Positive: Plvap, Emcn, Cdh5, Kdr, Egfl7, Mmrn2, Adgrf5, Esam, Pcdh17, Cyyr1
## Fabp4, Flt1, Tie1, Adgrl4, Myct1, Mcam, Ptprb, Prex2, Apold1, Gpihbp1
## Ramp2, Robo4, Epas1, Eng, Rasip1, Cd34, Fam167b, Grrp1, Arhgef15, Inhbb
## Negative: Dcn, Loxl1, Rpl32, Ccdc80, Serpinf1, Lox, Rpl21, Cpxm1, Col3a1, Col5a1
## Mfap5, Rps15a, Medag, Slit3, Pdgfra, Rps18, Rpl34, Col6a2, Lgi2, Cdh11
## Mmp2, Mmp23, Lgals1, Sfrp1, Fndc1, Rps2, Col5a2, Col6a1, Rcn3, Col1a2
## PC_ 5
## Positive: Ambp, Cldn3, Dcdc2a, Tstd1, Epcam, Cldn7, Pigr, Ces1d, Onecut1, Prox1
## Ttr, Cp, Fxyd3, Ehf, Slc28a3, 5330417C22Rik, Smim22, Tmem229a, Atp1b1, Hpn
## Pkhd1, Slc5a1, Fgfr3, Wfdc15b, Nr5a2, Bche, Gsta3, Fmo2, Gstt1, Scara3
## Negative: Tmsb10, Top2a, F2r, Ube2c, Gimap4, Pclaf, Inpp4b, Thy1, Ptpn22, Tpx2
## Ccna2, Trbc2, Rrm2, Has2, Birc5, Smc4, Cd3g, Cd3e, Tuba1a, Ms4a4b
## Bcl2l15, Cenpe, Cenpa, S1pr1, Kif11, Cenpf, Hist1h2ap, Racgap1, Dynap, Hmgb2
ElbowPlot(seurat_object)
seurat_object <- RunUMAP(seurat_object, reduction='pca', dims = 1:20)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 10:35:50 UMAP embedding parameters a = 0.9922 b = 1.112
## 10:35:50 Read 48154 rows and found 20 numeric columns
## 10:35:50 Using Annoy for neighbor search, n_neighbors = 30
## 10:35:50 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 10:35:54 Writing NN index file to temp file /var/folders/8b/2_66jk7x0jl1qc4177kp8svm0000gq/T//RtmpdAYuqi/file474d3f45c7ff
## 10:35:54 Searching Annoy index using 1 thread, search_k = 3000
## 10:36:03 Annoy recall = 100%
## 10:36:07 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 10:36:08 Initializing from normalized Laplacian + noise (using RSpectra)
## 10:36:10 Commencing optimization for 200 epochs, with 2070706 positive edges
## 10:36:22 Optimization finished
seurat_object <- FindNeighbors(seurat_object, reduction = 'pca', dims = 1:20)
## Computing nearest neighbor graph
## Computing SNN
seurat_object <- FindClusters(seurat_object, resolution = c(0.4))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 48154
## Number of edges: 1648610
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9476
## Number of communities: 23
## Elapsed time: 7 seconds
DimPlot(seurat_object, reduction='pca', label = TRUE)
DimPlot(seurat_object, label = TRUE)
Idents(seurat_object) <- 'RNA_snn_res.0.4'
p1 <- DimPlot(seurat_object, reduction = 'umap')
p2 <- DimPlot(seurat_object, reduction = 'umap', group.by = c("orig.ident"))
p1+p2
# after preprocessing (standard seurat workflow)
seurat_object <- IntegrateLayers(object = seurat_object, method = CCAIntegration, assay = "RNA", orig.reduction = "pca", new.reduction = "integrated.cca",
scale.layer = "scale.data",verbose = FALSE)
# follow up with dimensional reduction (again)
seurat_object <- RunUMAP(seurat_object, reduction='integrated.cca', dims = 1:20)
## 11:19:49 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:19:49 Read 48154 rows and found 20 numeric columns
## 11:19:49 Using Annoy for neighbor search, n_neighbors = 30
## 11:19:49 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:19:51 Writing NN index file to temp file /var/folders/8b/2_66jk7x0jl1qc4177kp8svm0000gq/T//RtmpdAYuqi/file474d7cc535dd
## 11:19:52 Searching Annoy index using 1 thread, search_k = 3000
## 11:20:01 Annoy recall = 100%
## 11:20:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:20:03 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:20:05 Commencing optimization for 200 epochs, with 2142582 positive edges
## 11:20:17 Optimization finished
seurat_object <- FindNeighbors(seurat_object, reduction = 'integrated.cca', dims = 1:20)
## Computing nearest neighbor graph
## Computing SNN
seurat_object <- FindClusters(seurat_object, resolution = c(0.6))
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 48154
## Number of edges: 1761701
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9197
## Number of communities: 23
## Elapsed time: 8 seconds
DimPlot(seurat_object, label = TRUE)
p1 <- DimPlot(seurat_object, reduction = 'umap', label = T)
p2 <- DimPlot(seurat_object, reduction = 'umap', group.by = c("orig.ident"))
p1+p2
# fibroblast + CAF markers
VlnPlot(seurat_object, features = c("Col1a1", "Col1a2", "Col3a1",
"Ngfr", "Pdgfra", "Acta2", "Cd74"))
FeaturePlot(seurat_object, features = c("Col1a1", "Col1a2", "Col3a1",
"Ngfr", "Pdgfra", "Acta2", "Cd74", "Plvap"))
# macrophage markers
VlnPlot(seurat_object, features = c("C1qa", "Lyz2",
"C1qb", "Mafb",
"C1qc", "Apoe",
"Csf1r", "Mrc1"))
FeaturePlot(seurat_object, features = c("C1qa", "Lyz2",
"C1qb", "Mafb",
"C1qc", "Apoe",
"Csf1r", "Mrc1"))
# other immune cells
# monocytes
VlnPlot(seurat_object, features = c("Ccr2", "Ly6c2"))
FeaturePlot(seurat_object, features = c("Ccr2", "Ly6c2"))
# cDCs
VlnPlot(seurat_object, features = c("Flt3", "Ccr7"))
FeaturePlot(seurat_object, features = c("Flt3", "Ccr7"))
# neutrophils
VlnPlot(seurat_object, features = c("Cxcr2", "Csf3r"))
FeaturePlot(seurat_object, features = c("Cxcr2", "Mpo", "Csf3r", "Elane"))
# pDCs
VlnPlot(seurat_object, features = c("Siglech", "Ccr9"))
FeaturePlot(seurat_object, features = c("Siglech", "Ccr9"))
# B cells
VlnPlot(seurat_object, features = c("Cd79b", "Jchain", "Mzb1"))
FeaturePlot(seurat_object, features = c("Cd79b", "Jchain", "Mzb1"))
# NK cells
VlnPlot(seurat_object, features = c("Gzma", "Prf1"))
FeaturePlot(seurat_object, features = c("Gzma", "Prf1"))
# T cells
VlnPlot(seurat_object, features = c("Lef1", "Il2rb", "Cd3e", "Gzmb", "Cd8a", "Xcl1"))
FeaturePlot(seurat_object, features = c("Lef1", "Il2rb", "Cd3e", "Gzmb", "Cd8a", "Xcl1"))
# tumor markers
VlnPlot(seurat_object, features = c("Epcam", "Krt8", "Krt7", "Krt19", "Top2a", "Gata6", "Hmga2"))
FeaturePlot(seurat_object, features = c("Epcam", "Krt8", "Krt7", "Krt19", "Top2a", "Hmga2", "Gata6"))
# acinar, ductal and ADM
VlnPlot(seurat_object, features = c("Krt19", "Cftr", "Cpb1", "Cpa1", "Cela2a", "Sox9", "Spp1", "Reg3a", "Crp"))
FeaturePlot(seurat_object, features = c("Krt19", "Cftr", "Cpb1", "Cpa1", "Cela2a", "Sox9", "Spp1", "Reg3a", "Crp"))
# endothelial cells
VlnPlot(seurat_object, features = c("Plvap", "Cdh5"))
FeaturePlot(seurat_object, features = c("Plvap", "Cdh5"))
# rename identified clusters
seurat_object <- RenameIdents(object = seurat_object, `10` = "Fibroblasts_1")
seurat_object <- RenameIdents(object = seurat_object, `14` = "Fibroblasts_2")
seurat_object <- RenameIdents(object = seurat_object, `20` = "Fibroblasts_3")
seurat_object <- RenameIdents(object = seurat_object, `1` = "Macrophages_1") # includes monocytes
seurat_object <- RenameIdents(object = seurat_object, `2` = "Macrophages_2")
seurat_object <- RenameIdents(object = seurat_object, `16` = "Macrophages_3")
seurat_object <- RenameIdents(object = seurat_object, `15` = "Ductal/ADM cells")
seurat_object <- RenameIdents(object = seurat_object, `13` = "Acinar cells")
seurat_object <- RenameIdents(object = seurat_object, `0` = "Tumor cells_1")
seurat_object <- RenameIdents(object = seurat_object, `4` = "Tumor cells_2")
seurat_object <- RenameIdents(object = seurat_object, `8` = "Tumor cells_3")
seurat_object <- RenameIdents(object = seurat_object, `9` = "Tumor cells_4")
seurat_object <- RenameIdents(object = seurat_object, `11` = "cDCs")
seurat_object <- RenameIdents(object = seurat_object, `21` = "pDCs")
seurat_object <- RenameIdents(object = seurat_object, `5` = "Neutrophils_1")
seurat_object <- RenameIdents(object = seurat_object, `18` = "Neutrophils_2")
seurat_object <- RenameIdents(object = seurat_object, `22` = "Neutrophils_3")
seurat_object <- RenameIdents(object = seurat_object, `3` = "Bcells_1")
seurat_object <- RenameIdents(object = seurat_object, `19` = "Bcells_2")
seurat_object <- RenameIdents(object = seurat_object, `6` = "NKcells")
seurat_object <- RenameIdents(object = seurat_object, `7` = "Tcells")
seurat_object <- RenameIdents(object = seurat_object, `12` = "Endothel")
DimPlot(seurat_object, reduction = "umap", label = T)
# save for easier import
saveRDS(seurat_object, file = "./GSE217846_RAW/data/timecourseKPC_preprocessed.rds")
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS 15.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Stockholm
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scDblFinder_1.18.0 SingleCellExperiment_1.26.0
## [3] SummarizedExperiment_1.34.0 Biobase_2.64.0
## [5] GenomicRanges_1.56.2 GenomeInfoDb_1.40.1
## [7] IRanges_2.38.1 S4Vectors_0.42.1
## [9] BiocGenerics_0.50.0 MatrixGenerics_1.16.0
## [11] matrixStats_1.5.0 lubridate_1.9.3
## [13] forcats_1.0.0 stringr_1.5.1
## [15] purrr_1.0.4 readr_2.1.5
## [17] tidyr_1.3.1 tibble_3.2.1
## [19] ggplot2_3.5.1 tidyverse_2.0.0
## [21] dplyr_1.1.4 Seurat_5.2.1
## [23] SeuratObject_5.0.2 sp_2.2-0
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.22 splines_4.4.1
## [3] later_1.4.1 BiocIO_1.14.0
## [5] bitops_1.0-9 polyclip_1.10-7
## [7] XML_3.99-0.17 fastDummies_1.7.5
## [9] lifecycle_1.0.4 edgeR_4.2.2
## [11] globals_0.16.3 lattice_0.22-6
## [13] MASS_7.3-61 magrittr_2.0.3
## [15] limma_3.60.6 plotly_4.10.4
## [17] sass_0.4.9 rmarkdown_2.29
## [19] jquerylib_0.1.4 yaml_2.3.10
## [21] metapod_1.12.0 httpuv_1.6.15
## [23] sctransform_0.4.1 spam_2.11-1
## [25] spatstat.sparse_3.1-0 reticulate_1.40.0
## [27] cowplot_1.1.3 pbapply_1.7-2
## [29] RColorBrewer_1.1-3 abind_1.4-8
## [31] zlibbioc_1.50.0 Rtsne_0.17
## [33] RCurl_1.98-1.16 GenomeInfoDbData_1.2.12
## [35] ggrepel_0.9.6 irlba_2.3.5.1
## [37] listenv_0.9.1 spatstat.utils_3.1-2
## [39] goftest_1.2-3 RSpectra_0.16-2
## [41] dqrng_0.4.1 spatstat.random_3.3-2
## [43] fitdistrplus_1.2-2 parallelly_1.42.0
## [45] DelayedMatrixStats_1.26.0 codetools_0.2-20
## [47] DelayedArray_0.30.1 scuttle_1.14.0
## [49] tidyselect_1.2.1 UCSC.utils_1.0.0
## [51] farver_2.1.2 viridis_0.6.5
## [53] ScaledMatrix_1.12.0 spatstat.explore_3.3-4
## [55] GenomicAlignments_1.40.0 jsonlite_1.9.0
## [57] BiocNeighbors_1.22.0 progressr_0.15.1
## [59] scater_1.32.1 ggridges_0.5.6
## [61] survival_3.7-0 tools_4.4.1
## [63] ica_1.0-3 Rcpp_1.0.14
## [65] glue_1.8.0 gridExtra_2.3
## [67] SparseArray_1.4.8 xfun_0.51
## [69] withr_3.0.2 fastmap_1.2.0
## [71] bluster_1.14.0 digest_0.6.37
## [73] rsvd_1.0.5 timechange_0.3.0
## [75] R6_2.6.1 mime_0.12
## [77] colorspace_2.1-1 scattermore_1.2
## [79] tensor_1.5 spatstat.data_3.1-4
## [81] generics_0.1.3 data.table_1.16.4
## [83] rtracklayer_1.64.0 httr_1.4.7
## [85] htmlwidgets_1.6.4 S4Arrays_1.4.1
## [87] uwot_0.2.2 pkgconfig_2.0.3
## [89] gtable_0.3.6 lmtest_0.9-40
## [91] XVector_0.44.0 htmltools_0.5.8.1
## [93] dotCall64_1.2 scales_1.3.0
## [95] png_0.1-8 spatstat.univar_3.1-1
## [97] scran_1.32.0 knitr_1.49
## [99] rstudioapi_0.17.1 tzdb_0.4.0
## [101] reshape2_1.4.4 rjson_0.2.23
## [103] nlme_3.1-166 curl_6.2.1
## [105] cachem_1.1.0 zoo_1.8-12
## [107] KernSmooth_2.23-24 vipor_0.4.7
## [109] parallel_4.4.1 miniUI_0.1.1.1
## [111] ggrastr_1.0.2 restfulr_0.0.15
## [113] pillar_1.10.1 grid_4.4.1
## [115] vctrs_0.6.5 RANN_2.6.2
## [117] promises_1.3.2 BiocSingular_1.20.0
## [119] beachmat_2.20.0 xtable_1.8-4
## [121] cluster_2.1.6 beeswarm_0.4.0
## [123] evaluate_1.0.3 locfit_1.5-9.11
## [125] cli_3.6.4 compiler_4.4.1
## [127] Rsamtools_2.20.0 rlang_1.1.5
## [129] crayon_1.5.3 future.apply_1.11.3
## [131] labeling_0.4.3 ggbeeswarm_0.7.2
## [133] plyr_1.8.9 stringi_1.8.4
## [135] viridisLite_0.4.2 deldir_2.0-4
## [137] BiocParallel_1.38.0 munsell_0.5.1
## [139] Biostrings_2.72.1 lazyeval_0.2.2
## [141] spatstat.geom_3.3-5 Matrix_1.7-1
## [143] RcppHNSW_0.6.0 hms_1.1.3
## [145] patchwork_1.3.0 sparseMatrixStats_1.16.0
## [147] future_1.34.0 statmod_1.5.0
## [149] shiny_1.10.0 ROCR_1.0-11
## [151] igraph_2.1.4 bslib_0.9.0
## [153] xgboost_1.7.8.1